
#include "gsl_ext_function.h"
#include "gsl_ext_testfunctions.h"

#include <math.h>
#include <string.h>

#define INCREASE_EVAL_COUNTER(p) ((gsl_ext_multimin_function_params *)p)->f_eval_counter++

#define GET_M ((gsl_ext_testfunction_params *)((gsl_ext_multimin_function_params *)params)->params)->m
#define GET_N ((gsl_ext_testfunction_params *)((gsl_ext_multimin_function_params *)params)->params)->n

#define EVAL_COUNTING

static double fx;
static int i;

static const double gaussian_y[] = { 0.0009, 0.0044, 0.0175, 0.0540, 0.1295, 0.2420, 0.3521, 0.3989, 
                                     0.3521, 0.2420, 0.1295, 0.0540, 0.0175, 0.0044, 0.0009 };

static int testfunction_params_no_check(const gsl_ext_testfunction_params *p)
{
  return 1;
}

int biggsexp6_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->m >= p->n);
}

int browndennis_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->m >= p->n);
}

int chebyquad_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->m >= p->n);
}

int extendedpowellsingular_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->n % 4 == 0);
}

int extendedrosenbrock_check_params(const gsl_ext_testfunction_params *p)
{
  if(p->n % 2 != 0)
    return 0;
	
  return 1;
}

int gulf_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->m >= p->n && p->m <= 100);
}

int watson_check_params(const gsl_ext_testfunction_params *p)
{
  return (p->n >= 2 && p->n <= 31);
}

int wood_check_params(const gsl_ext_testfunction_params *p)
{
  return 1;
}

double beale_f(const gsl_vector * x, void *params)
{
  fx = 0.0;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
	
  double f1 = (1.5 - x1 * (1.0 - x2));
  double f2 = (2.25 - x1 * (1.0 - x2*x2));
  double f3 = (2.625 - x1 * (1.0 - x2*x2*x2));
	
  fx += f1 * f1;
  fx += f2 * f2;
  fx += f3 * f3;
   
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double biggsexp6_f(const gsl_vector *x, void *params)
{
  fx = 0.0;
  double fix;
  double ti;
  double yi;
	
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
  double x4 = gsl_vector_get(x, 3);
  double x5 = gsl_vector_get(x, 4);
  double x6 = gsl_vector_get(x, 5);
	
  for(i = 1; i <= GET_M; i++)
  {
    ti = 0.1 * i;
    yi = exp(-ti) - 5.0*exp(-10.0*ti) + 3.0*exp(-4.0*ti);
    fix = x3*exp(-ti*x1) - x4*exp(-ti*x2) + x6*exp(-ti*x5) - yi;
    fx += fix * fix;
  }
	
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double box_f(const gsl_vector * x, void *params)
{
  double fix;
  double ti;
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
	
  fx = 0.0;
  
  for(i = 1; i <= GET_M; i++)
  {
    ti = 0.1 * i;
    fix = exp(-ti * x1) - exp(-ti * x2) - x3 * (exp(-ti) - exp(-10.0 * ti));
    fx += fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double brownbadlyscaled_f(const gsl_vector *x, void *params)
{
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  
  double f1 = x1 - 1e6;
  double f2 = x2 - 2.0 * 1e-6;
  double f3 = x1 * x2 - 2.0;
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return f1 * f1 + f2 * f2 + f3 * f3;
}

double browndennis_f(const gsl_vector *x, void *params)
{
  double ti;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
  double x4 = gsl_vector_get(x, 3);
  
  double fix;
  double term1, term2;
  
  fx = 0.0;
  
  for(i = 1; i <= GET_M; i++)
  {
    ti = i / 5.0;
    
    term1 = x1 + ti * x2 - exp(ti);
    term2 = x3 + x4 * sin(ti) - cos(ti);
    
    fix = term1 * term1 + term2 * term2;
    fx += fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double chebyquad_f(const gsl_vector *x, void *params)
{
  double fix;
  int i, j, k;
  double sum;
  double Ti, Tim1, prevTi;
  
  fx = 0.0;
	
  for(i = 1; i <= GET_M; i++)
  {
    sum = 0.0;
    for(j = 1; j <= GET_N; j++)
    {
      Tim1 = 1.0;
      Ti = gsl_vector_get(x, j - 1);
      for(k = 2; k <= i; k++)
      {
        prevTi = Ti;
        Ti = 2.0 * gsl_vector_get(x, j - 1) * Ti - Tim1;
        Tim1 = prevTi;
      }
      
      sum += Ti;
    }
    sum /= GET_N;
  
    fix = sum;
  
    if(i % 2 == 0)
      fix += 1.0 / (i * i - 1.0);
  
    fx += fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double extendedpowellsingular_f(const gsl_vector * x, void *params)
{
  double fix;
  double x1, x2, x3, x4;
  
  fx = 0.0;
  
  for(i = 0; i < GET_N; i += 4)
  {
    x1 = gsl_vector_get(x, i);
    x2 = gsl_vector_get(x, i + 1);
    x3 = gsl_vector_get(x, i + 2);
    x4 = gsl_vector_get(x, i + 3);
	      
    fix = x1 + 10.0 * x2;
    fx += fix * fix;
  
    fix = x3 - x4;
    fx += 5.0 * fix * fix;
    
    fix = (x2 - 2.0 * x3) * (x2 - 2.0 * x3);
    fx += fix * fix;
  
    fix = (x1 - x4) * (x1 - x4);
    fx += 10.0 * fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double extendedrosenbrock_f(const gsl_vector * x, void *params)
{
  fx = 0.0;
  
  double f1, f2;
  double x1, x2;
  
  for(i = 0; i < GET_N; i += 2)
  {
    x1 = gsl_vector_get(x, i);
    x2 = gsl_vector_get(x, i + 1);
    
    f1 = 10.0 * (x2 - x1*x1);
    f2 = 1.0 - x1;
    
    fx += f1*f1 + f2*f2;
  }
	
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double gaussian_f(const gsl_vector * x, void *params)
{
  double fix;
  double ti;
  
  fx = 0.0;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
  
  for(i = 1; i <= 15; i++)
  {
    ti = (8.0 - i) / 2.0;
    fix = x1 * exp(-x2 * (ti-x3) * (ti-x3) / 2.0) - gaussian_y[i - 1];
    fx += fix * fix;
  }

  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double gulf_f(const gsl_vector * x, void *params)
{
  double fix;
  double ti;
  double yi;
  
  fx = 0.0;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
	
  for(i = 1; i <= GET_M; i++)
  {
    ti = i / 100.0;
    yi = 25.0 + pow(-50.0 * log(ti), 2.0 / 3.0);
    fix = exp(-pow(fabs(yi - x2), x3) / x1) - ti;
    fx += fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double helicalvalley_f(const gsl_vector * x, void *params)
{
  double phi;
        
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
	
  double second_term;
	
  if(x1 > 0.0)
    phi = 1.0 / (2.0 * M_PI) * atan(x2 / x1);
  else
    phi = 1.0 / (2.0 * M_PI) * atan(x2 / x1) + 0.5;
	
  second_term = 10.0 * (sqrt(x1*x1 + x2*x2) - 1.0);
	
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return 100.0 * (x3 - 10.0*phi) * (x3 - 10.0*phi) + 
         second_term * second_term + 
         x3 * x3;
}

double penaltyfunctionI_f(const gsl_vector * x, void *params)
{
  fx = 0.0;
	
  double fx = 0.0;
  double fix;
  int i;
  double sum = 0.0;
  
  for(i = 1; i <= GET_N; i++)
  {
    fix = gsl_vector_get(x, i - 1) - 1.0;
    fx += 1e-5 * fix * fix;
  }
  
  for(i = 1; i <= GET_N; i++)
    sum += gsl_vector_get(x, i - 1) * gsl_vector_get(x, i - 1);
  
  fx += (sum - 0.25) * (sum - 0.25);
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double penaltyfunctionII_f(const gsl_vector * x, void *params)
{
  double fix;
  double sum = 0.0;
  double yi;
  
  fx = 0.0;
  
  fx += (gsl_vector_get(x, 0) - 0.2) * (gsl_vector_get(x, 0) - 0.2);
  
  for(i = 2; i < 2*GET_N; i++)
  {
    yi = exp(i / 10.0) + exp((i - 1) / 10.0);
  
    if(i <= GET_N)
      fix = exp(gsl_vector_get(x, i - 1) / 10.0) + exp(gsl_vector_get(x, i - 2) / 10.0) - yi;
    else
      fix = exp(gsl_vector_get(x, i - GET_N) / 10.0) - exp(-1.0 / 10.0);
  
    fx += 1e-5 * fix * fix;
  }
  
  for(i = 1; i <= GET_N; i++)
    sum += (GET_N - i + 1) * gsl_vector_get(x, i - 1) * gsl_vector_get(x, i - 1);
	
  fx += (sum - 1.0) * (sum - 1.0);
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double powellbadlyscaled_f(const gsl_vector * x, void *params)
{
  fx = 0.0;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
	
  double fx1 = 10000.0 * x1 * x2 - 1.0;
  double fx2 = exp(-x1) + exp(-x2) - 1.0001;
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx1 * fx1 + fx2 * fx2;
}

double trigonometric_f(const gsl_vector * x, void *params)
{
  double fix;
  int i, j;
  
  fx = 0.0;
  
  for(i = 1; i <= GET_N; i++)
  {
    fix = 0.0;
	      
    fix += GET_N;
    for(j = 1; j <= GET_N; j++)
      fix -= cos(gsl_vector_get(x, j - 1));
    fix += i * (1.0 - cos(gsl_vector_get(x, i - 1))) - sin(gsl_vector_get(x, i - 1));
          
    fx += fix * fix;
  }
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double variablydimensioned_f(const gsl_vector * x, void *params)
{
  double sum = 0.0;
  
  fx = 0.0;
  
  for(i = 1; i <= GET_N; i++)
    fx += (gsl_vector_get(x, i - 1) - 1.0) * (gsl_vector_get(x, i - 1) - 1.0);
  
  for(i = 1; i <= GET_N; i++)
    sum += i * (gsl_vector_get(x, i - 1) - 1.0);
  
  fx += sum * sum;
  
  fx += sum * sum * sum * sum;
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double watson_f(const gsl_vector * x, void *params)
{
  double fix;
  int i, j;
  double s;
  double ti;
  
  fx = 0.0;
  
  for(i = 1; i <= 29; i++)
  {
    fix = 0.0;
    
    ti = ((double)i) / 29.0;
    s = 0.0;
    for(j = 2; j <= GET_N; j++)
      s += (j - 1) * gsl_vector_get(x, j - 1) * pow(ti, j - 2);
    fix += s;
    
    s = 0.0;
    for(j = 1; j <= GET_N; j++)
      s += gsl_vector_get(x, j - 1) * pow(ti, j - 1);
    fix -= s *s + 1.0;
    
    fx += fix * fix;
  }
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  
  fx += x1 * x1;
  fx += (x2 - x1 * x1 - 1.0) * (x2 - x1 * x1 - 1.0);
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

double wood_f(const gsl_vector * x, void *params)
{
  fx = 0.0;
  
  double x1 = gsl_vector_get(x, 0);
  double x2 = gsl_vector_get(x, 1);
  double x3 = gsl_vector_get(x, 2);
  double x4 = gsl_vector_get(x, 3);
  
  fx += 100.0 * (x2 - x1 * x1) * (x2 - x1 * x1);
  fx += (1.0 - x1) * (1.0 - x1);
  fx += 90.0 * (x4 - x3 * x3) * (x4 - x3 * x3);
  fx += (1.0 - x3) * (1.0 - x3);
  fx += 10.0 * (x2 + x4 - 2.0) * (x2 + x4 - 2.0);
  fx += 0.1 * (x2 - x4) * (x2 - x4);
  
  #ifdef EVAL_COUNTING
  INCREASE_EVAL_COUNTER(params);
  #endif
  
  return fx;
}

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_beale = 
{
  "Beale", 
  &beale_f, 
  2,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_biggsexp6 = 
{
  "Biggs EXP6", 
  &biggsexp6_f, 
  6,
  biggsexp6_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_box = 
{
  "Box", 
  &box_f, 
  3, 
  biggsexp6_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_brownbadlyscaled = 
{
  "Brown badly scaled", 
  &brownbadlyscaled_f, 
  2,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_browndennis = 
{
  "Brown and Dennis", 
  &browndennis_f, 
  4,
  browndennis_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_chebyquad = 
{
  "Chebyquad", 
  &chebyquad_f, 
  0,
  chebyquad_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_extendedpowellsingular = 
{
  "Extended Powell singular", 
  &extendedpowellsingular_f, 
  0,
  extendedpowellsingular_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_extendedrosenbrock = 
{
  "Extended Rosenbrock", 
  &extendedrosenbrock_f, 
  0,
  extendedrosenbrock_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_gaussian = 
{
  "Gaussian", 
  &gaussian_f, 
  3,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_gulf = 
{
  "Gulf", 
  &gulf_f, 
  3,
  gulf_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_helicalvalley = 
{
  "Gaussian", 
  &helicalvalley_f, 
  3,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_penaltyfunctionI = 
{
  "Penalty function I", 
  &penaltyfunctionI_f, 
  0,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_penaltyfunctionII = 
{
  "Penalty function II", 
  &penaltyfunctionII_f, 
  0,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_powellbadlyscaled = 
{
  "Powell badly scaled", 
  &powellbadlyscaled_f, 
  2,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_trigonometric = 
{
  "Trigonometric", 
  &trigonometric_f, 
  0,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_variablydimensioned = 
{
  "Variably dimensioned", 
  &variablydimensioned_f, 
  0,
  testfunction_params_no_check
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_watson = 
{
  "Watson", 
  &watson_f, 
  0,
  watson_check_params
};

static const gsl_ext_testfunction_type gsl_ext_testfunction_type_wood = 
{
  "Wood", 
  &wood_f, 
  4,
  wood_check_params
};

const gsl_ext_testfunction_type *
gsl_ext_testfunction_beale = &gsl_ext_testfunction_type_beale;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_biggsexp6 = &gsl_ext_testfunction_type_biggsexp6;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_box = &gsl_ext_testfunction_type_box;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_brownbadlyscaled = &gsl_ext_testfunction_type_brownbadlyscaled;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_browndennis = &gsl_ext_testfunction_type_browndennis;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_chebyquad = &gsl_ext_testfunction_type_chebyquad;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_extendedpowellsingular = &gsl_ext_testfunction_type_extendedpowellsingular;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_extendedrosenbrock = &gsl_ext_testfunction_type_extendedrosenbrock;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_gaussian = &gsl_ext_testfunction_type_gaussian;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_gulf = &gsl_ext_testfunction_type_gulf;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_helicalvalley = &gsl_ext_testfunction_type_helicalvalley;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_penaltyfunctionI = &gsl_ext_testfunction_type_penaltyfunctionI;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_penaltyfunctionII = &gsl_ext_testfunction_type_penaltyfunctionII;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_powellbadlyscaled = &gsl_ext_testfunction_type_powellbadlyscaled;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_trigonometric = &gsl_ext_testfunction_type_trigonometric;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_variablydimensioned = &gsl_ext_testfunction_type_variablydimensioned;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_watson = &gsl_ext_testfunction_type_watson;
const gsl_ext_testfunction_type *
gsl_ext_testfunction_wood = &gsl_ext_testfunction_type_wood;

int gsl_ext_testfunction_check_params(const gsl_ext_testfunction_type *T, const gsl_ext_testfunction_params *p)
{
  return T->check_params(p);
}

int gsl_ext_testfunction_needs_params(const gsl_ext_testfunction_type *T)
{
  return (T->check_params != &testfunction_params_no_check);
}

const gsl_ext_testfunction_type *gsl_ext_get_testfunction(char *expr)
{
  if(strcmp(expr, "beale") == 0)
    return gsl_ext_testfunction_beale;
  else if(strcmp(expr, "biggsexp6") == 0)
    return gsl_ext_testfunction_biggsexp6;
  else if(strcmp(expr, "box") == 0)
    return gsl_ext_testfunction_box;
  else if(strcmp(expr, "brownbadlyscaled") == 0)
    return gsl_ext_testfunction_brownbadlyscaled;
  else if(strcmp(expr, "browndennis") == 0)
    return gsl_ext_testfunction_browndennis;
  else if(strcmp(expr, "chebyquad") == 0)
    return gsl_ext_testfunction_chebyquad;
  else if(strcmp(expr, "extendedpowellsingular") == 0)
    return gsl_ext_testfunction_extendedpowellsingular;
  else if(strcmp(expr, "extendedrosenbrock") == 0)
    return gsl_ext_testfunction_extendedrosenbrock;
  else if(strcmp(expr, "gaussian") == 0)
    return gsl_ext_testfunction_gaussian;
  else if(strcmp(expr, "gulf") == 0)
    return gsl_ext_testfunction_gulf;
  else if(strcmp(expr, "helicalvalley") == 0)
    return gsl_ext_testfunction_helicalvalley;
  else if(strcmp(expr, "penaltyfunctionI") == 0)
    return gsl_ext_testfunction_penaltyfunctionI;
  else if(strcmp(expr, "penaltyfunctionII") == 0)
    return gsl_ext_testfunction_penaltyfunctionII;
  else if(strcmp(expr, "powellbadlyscaled") == 0)
    return gsl_ext_testfunction_powellbadlyscaled;
  else if(strcmp(expr, "trigonometric") == 0)
    return gsl_ext_testfunction_trigonometric;
  else if(strcmp(expr, "variablydimensioned") == 0)
    return gsl_ext_testfunction_variablydimensioned;
  else if(strcmp(expr, "watson") == 0)
    return gsl_ext_testfunction_watson;
  else if(strcmp(expr, "wood") == 0)
    return gsl_ext_testfunction_wood;
  else
    return NULL;
}
